In [ ]:
import geoengine as ge
from geoengine.workflow_builder.operators import GdalSource, TemporalRasterAggregation, RasterStacker, RenameBands, \
    Expression,RasterTypeConversion, RasterBandDescriptor
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from typing import Tuple, Optional
In [ ]:
ge.initialize("https://zentrale.app.geoengine.io/api", token="51c46772-6c55-4867-8e42-c9fd8b614820")
In [ ]:
def _query_rectangle(*,
                     center: Tuple[float, float],
                     time: np.datetime64,
                     radius_px: float = 512) -> ge.QueryRectangle:
    resolution = ge.SpatialResolution(10, 10)
    bbox = ge.BoundingBox2D(
        xmin=center[0] - resolution.x_resolution * radius_px,
        xmax=center[0] + resolution.x_resolution * radius_px,
        ymin=center[1] - resolution.y_resolution * radius_px,
        ymax=center[1] + resolution.y_resolution * radius_px,
    )
    return ge.QueryRectangle(
        spatial_bounds=bbox,
        time_interval=ge.TimeInterval(time),
        srs='EPSG:32632',
        resolution=resolution,
    )

async def compute_rgb_array(*,
                      center: Tuple[float, float],
                      time: np.datetime64,
                      radius_px: float = 512,
                      ) -> xr.DataArray:
    '''Compute RGB of area'''

    def cloud_free(band: GdalSource):
        return Expression(
            expression = 'if B > 2 && B < 8 { A } else { NODATA }',
            source=RasterStacker(
            sources=[
                band,
                RasterTypeConversion(
                    GdalSource("_:5779494c-f3a2-48b3-8a2d-5fbba8c5b6c5:`UTM32N:SCL`"),
                    output_data_type='U16',
                ),
            ],
            )
        )
    
    query_rectangle = _query_rectangle(
        center=center,
        time=time,
        radius_px=radius_px,
    )

    workflow = TemporalRasterAggregation(
        aggregation_type='mean', # 'percentileEstimate',
        # percentile=0.5,
        granularity='months',
        window_size=1,
        ignore_no_data=True,
        source=RasterStacker(
            sources=[
            cloud_free(GdalSource("_:5779494c-f3a2-48b3-8a2d-5fbba8c5b6c5:`UTM32N:B04`")),
            cloud_free(GdalSource("_:5779494c-f3a2-48b3-8a2d-5fbba8c5b6c5:`UTM32N:B03`")),
            cloud_free(GdalSource("_:5779494c-f3a2-48b3-8a2d-5fbba8c5b6c5:`UTM32N:B02`")),
            ],
            rename=RenameBands.rename(['red', 'green', 'blue']),
        )
    )

    workflow = ge.register_workflow(workflow)

    data_array = await workflow.raster_stream_into_xarray(
        query_rectangle=query_rectangle,
        clip_to_query_rectangle=True,
        bands=[0,1,2], # TODO: improve for user, default = all? where are the band names?
    )

    return data_array


async def compute_water_bodies(*,
                      center: Tuple[float, float],
                      time: np.datetime64,
                      radius_px: float = 512,
                      water_body_time_aggregate_months: int = 1,
                      ) -> xr.DataArray:
    '''Compute water body classification'''
    query_rectangle = _query_rectangle(
        center=center,
        time=time,
        radius_px=radius_px,
    )

    workflow = TemporalRasterAggregation(
        aggregation_type='max',
        granularity='months',
        window_size=water_body_time_aggregate_months,
        ignore_no_data=True,
        source=Expression(
            expression='''
                let ndwi = (A-B)/(A+B);
                if C <= 2 || C >= 8 {
                    NODATA
                } else if ndwi >= 0.2 {
                    3
                } else if ndwi >= 0 {
                    2
                } else if ndwi >= -0.3 {
                    1
                } else {
                    0
                }
            ''',
            output_type='U8',
            # output_band=RasterBandDescriptor( # TODO: this is broken in the backend
            #     "water bodies",
            #     measurement=ge.ClassificationMeasurement(
            #         "water bodies",
            #         classes={
            #             0: 'Drought, non-aqueous surfaces',
            #             1: 'Moderate drought, non-aqueous surfaces',
            #             2: 'Flooding, humidity',
            #             3: 'Water surface',
            #         }
            #     ),
            # ),
            source=
                RasterStacker(
                    sources=[
                        GdalSource("_:5779494c-f3a2-48b3-8a2d-5fbba8c5b6c5:`UTM32N:B03`"),
                        GdalSource("_:5779494c-f3a2-48b3-8a2d-5fbba8c5b6c5:`UTM32N:B08`"),
                        RasterTypeConversion(
                            GdalSource("_:5779494c-f3a2-48b3-8a2d-5fbba8c5b6c5:`UTM32N:SCL`"),
                            output_data_type='U16',
                        ),
                    ],
                ),
        ),
    )

    workflow = ge.register_workflow(workflow)

    data_array = await workflow.raster_stream_into_xarray(
        query_rectangle=query_rectangle,
        clip_to_query_rectangle=True,
        bands=[0], # TODO: improve for user, default = all? where are the band names?
    )

    return data_array


def plot(rbg_array: xr.DataArray, water_bodies_array: Optional[xr.DataArray] = None, height_px: int = 1024):
    px = 1/plt.rcParams['figure.dpi']  # pixel in inches

    fig, axes = plt.subplot_mosaic("AB")

    fig.set_figheight(height_px*px)
    fig.set_figwidth(2 * height_px*px)

    rbg_array.isel(time=0).plot.imshow(
        rgb="band",
        vmax=4000,
        ax=axes['A'],
    )

    if water_bodies_array is None:
        fig.show()
        return
    
    water_bodies_array.isel(band=0, time=0, drop=True).plot.imshow(
        colors=[plt.cm.Paired.colors[i] for i in [2, 3, 0, 1]],
        levels=5,
        vmin=0,
        vmax=4,
        ax=axes['B']
    )

    axes['B'].images[-1].colorbar.set_ticks(
        [0,1,2,3],
        labels=[
            'Drought, non-aqueous surfaces',
            'Moderate drought, non-aqueous surfaces',
            'Flooding, humidity',
            'Water surface',
        ],
    )

    fig.show()

    fig, axes = plt.subplot_mosaic("C")

    fig.set_figheight(height_px*px)
    fig.set_figwidth(height_px*px)

    rbg_array.isel(time=0).plot.imshow(
        rgb="band",
        vmax=4000,
        ax=axes['C'],
    )

    water_bodies_array.isel(band=0, time=0, drop=True).plot.imshow(
        colors=[(0,0,0,0), plt.cm.Paired.colors[1]],
        levels=2,
        vmin=3,
        vmax=4,
        add_colorbar=False,
        ax=axes['C']
    )

    fig.show()
In [ ]:
marburg_center_utm = [483843, 5627614]
In [ ]:
marburg_rgb = await compute_rgb_array(center=marburg_center_utm, time=np.datetime64("2022-07-01T00:00:00"))
/home/beilschmidt/git/geoengine-python/env/lib/python3.10/site-packages/rasterio/windows.py:314: RasterioDeprecationWarning: The height, width, and precision parameters are unused, deprecated, and will be removed in 2.0.0.
  warnings.warn(
In [ ]:
marburg_water = await compute_water_bodies(center=marburg_center_utm, time=np.datetime64("2022-07-01T00:00:00"))
/home/beilschmidt/git/geoengine-python/env/lib/python3.10/site-packages/rasterio/windows.py:314: RasterioDeprecationWarning: The height, width, and precision parameters are unused, deprecated, and will be removed in 2.0.0.
  warnings.warn(
In [ ]:
plot(rbg_array=marburg_rgb, water_bodies_array=marburg_water)
/tmp/ipykernel_9616/3560293408.py:180: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure.
  fig.show()
/tmp/ipykernel_9616/3560293408.py:202: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure.
  fig.show()
/home/beilschmidt/git/geoengine-python/env/lib/python3.10/site-packages/matplotlib/cm.py:478: RuntimeWarning: invalid value encountered in cast
  xx = (xx * 255).astype(np.uint8)
No description has been provided for this image
No description has been provided for this image
In [ ]:
weimar_center_utm = [479763, 5623423]

weimar_rgb = await compute_rgb_array(center=weimar_center_utm, time=np.datetime64("2022-07-01T00:00:00"))

weimar_water = await compute_water_bodies(center=weimar_center_utm, time=np.datetime64("2022-07-01T00:00:00"))

plot(rbg_array=weimar_rgb, water_bodies_array=weimar_water)
/home/beilschmidt/git/geoengine-python/env/lib/python3.10/site-packages/rasterio/windows.py:314: RasterioDeprecationWarning: The height, width, and precision parameters are unused, deprecated, and will be removed in 2.0.0.
  warnings.warn(
/home/beilschmidt/git/geoengine-python/env/lib/python3.10/site-packages/rasterio/windows.py:314: RasterioDeprecationWarning: The height, width, and precision parameters are unused, deprecated, and will be removed in 2.0.0.
  warnings.warn(
/tmp/ipykernel_9616/3560293408.py:180: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure.
  fig.show()
/tmp/ipykernel_9616/3560293408.py:202: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure.
  fig.show()
/home/beilschmidt/git/geoengine-python/env/lib/python3.10/site-packages/matplotlib/cm.py:478: RuntimeWarning: invalid value encountered in cast
  xx = (xx * 255).astype(np.uint8)
No description has been provided for this image
No description has been provided for this image
In [ ]:
koeln_center_utm = [356766, 5644819]

koeln_rgb = await compute_rgb_array(center=koeln_center_utm, time=np.datetime64("2022-07-01T00:00:00"))

koeln_water = await compute_water_bodies(center=koeln_center_utm, time=np.datetime64("2022-07-01T00:00:00"))

plot(rbg_array=koeln_rgb, water_bodies_array=koeln_water)
/home/beilschmidt/git/geoengine-python/env/lib/python3.10/site-packages/rasterio/windows.py:314: RasterioDeprecationWarning: The height, width, and precision parameters are unused, deprecated, and will be removed in 2.0.0.
  warnings.warn(
/home/beilschmidt/git/geoengine-python/env/lib/python3.10/site-packages/rasterio/windows.py:314: RasterioDeprecationWarning: The height, width, and precision parameters are unused, deprecated, and will be removed in 2.0.0.
  warnings.warn(
/tmp/ipykernel_9616/3560293408.py:180: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure.
  fig.show()
/tmp/ipykernel_9616/3560293408.py:202: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure.
  fig.show()
No description has been provided for this image
No description has been provided for this image